Read in, format, and visualize data.

# Read in target data

target_data <- read.csv("terrestrial_30min-targets.csv")

# Read in gap-filled temperature data

temp_data <- read.csv("temp_gapfill.csv")

# Convert time column from character to POSIXct

# Time zone in target data is UTC according to https://docs.google.com/document/d/1l7sxBk-z-GHTlk50rdxP0lPTwJzFJ2gykclINkMsWcc/edit

target_data$time <- ymd_hms(target_data$time, tz = "UTC")
temp_data$startDateTime <- ymd_hms(temp_data$startDateTime, tz = "UTC")  # Ask Marcus to verify UTC

# Convert siteID column from character to factor

target_data$siteID <- as.factor(target_data$siteID)
temp_data$siteID <- as.factor(temp_data$siteID)

# Plot all data

ggplot(data = target_data) +
  geom_point(aes(x = time, y = nee, color = siteID)) +
  facet_grid(siteID ~ .) + 
  labs(title = "Full Time Period: 01 Feb 2017 - 31 Jan 2021",
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw() +
  theme(legend.position = "none")

# Plot data for example week (2020-06-07 through 2020-06-13)

ggplot(data = filter(target_data,
                     time >= "2020-06-07",
                     time < "2020-06-14")) +
  geom_point(aes(x = time, y = nee, color = siteID)) +
  facet_grid(siteID ~ .) + 
  labs(title = "Example Time Period: 07 Jun 2020 - 13 Jun 2020",
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw() +
  theme(legend.position = "none")

Construct null model (random walk) for NIMBLE.

RandomWalk <- nimbleCode({
  
  # Note:
  # x = real NEE (state variable)
  # y = observed NEE
  
  # Priors
  
  x[1] ~ dnorm(x_ic, sd = sd_ic)
  sd_add ~ dunif(0, 100)
  
  # Process model
  
  for(t in 2:n){
    pred[t] <- x[t-1]
    x[t] ~ dnorm(pred[t], sd = sd_add)
  }
  
  # Data model
  
  for(t in 1:n){
    y[t] ~ dnorm(x[t], sd = sd_obs[t])
  }
  
})

Construct temperature dynamic linear model (DLM) for NIMBLE.

TempDLM <- nimbleCode({
  
  # Note:
  # x = real NEE (state variable)
  # y = observed NEE
  
  # Priors
  
  x[1] ~ dnorm(x_ic, sd = sd_ic)
  sd_add ~ dunif(0, 100)
  b1 ~ dnorm(0, 100)  # Temperature effect coefficient
  
  # Process model
  
  for(t in 2:n){
    pred[t] <- x[t-1] + b1 * temp[t]
    x[t] ~ dnorm(pred[t], sd = sd_add)
  }
  
  # Data model
  
  for(t in 1:n){
    y[t] ~ dnorm(x[t], sd = sd_obs[t])
  }
  
})

Because of the size of the data set, a subset of the data had to be used for model fitting. The variables controlling the size of the data subset are specified below.

# Note: B/c of size of data set (and potentially number of NAs), need to use subset of data

subset_length_days <- 7  # Length of subset [d]

subset_length_timesteps <- subset_length_days * 24 * 2  # Number of time steps in subset (days x hours/day x half hours/hour)

Specify data from BART to be used for model fitting. (Because of the large difference between NEE values in the winter versus other times of year, the additional restriction of limiting the data subset to April-November was used.)

# Separate data by site

target_data_BART <- filter(target_data, siteID == "BART")
temp_data_BART <- filter(temp_data, siteID == "BART")

# Find data subset w/ fewest number of NAs

min_n_NAs <- subset_length_timesteps  # Initialize scalar for fewest NAs

BART_start_i <- NA  # Initialize BART starting index

# Loop through data to find data subset (not during winter) with minimum number of NAs for specified subset length

for(i in 1:(length(target_data_BART$nee) - subset_length_timesteps + 1)){
  
  num_NAs <- sum(is.na(target_data_BART$nee[i:(i + subset_length_timesteps - 1)]))  # Count number of NEE NAs in data subset

  if(num_NAs <= min_n_NAs & sum(substr(target_data_BART$time[i], 6, 7) != c("12", "01", "02", "03")) == 4){  # Note: Using <= so that most recent data is used; also, excluding data from December through March
    
    BART_start_i <- i  # Save new starting index
    min_n_NAs <- num_NAs  # Save new lowest number of NAs
    
  }
  
}

# Calculate data subset end index

BART_end_i <- BART_start_i + subset_length_timesteps - 1

# Subset data

target_data_BART_sub <- target_data_BART[BART_start_i:BART_end_i,]

# Interpolate NEE values to estimate sd_obs

nee_interp_BART_sub <- approx(x = target_data_BART_sub$time[!is.na(target_data_BART_sub$nee)],
                              y = target_data_BART_sub$nee[!is.na(target_data_BART_sub$nee)],
                              xout = target_data_BART_sub$time,
                              method = "linear",
                              yleft = mean(target_data_BART_sub$nee, na.rm = TRUE),
                              yright = mean(target_data_BART_sub$nee, na.rm = TRUE),
                              na.rm = TRUE)

sd_obs_BART_sub <- rep(NA, length(target_data_BART_sub$nee))  # Initialize sd_obs vector

for(i in 1:length(sd_obs_BART_sub)){
  
  if(nee_interp_BART_sub$y[i] >= 0){  # Calculate sd_obs if NEE is positive
    
    sd_obs_BART_sub[i] <- target_data_BART_sub$nee_sd_intercept[1] + target_data_BART_sub$ nee_sd_slopeP[1] * nee_interp_BART_sub$y[i]  # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
    
  }else{  # Calculate sd_obs if NEE is negative
    
    sd_obs_BART_sub[i] <- target_data_BART_sub$nee_sd_intercept[1] + target_data_BART_sub$ nee_sd_slopeN[1] * nee_interp_BART_sub$y[i]
    
  }
  
}

# Match temperature data to subset NEE data

BART_temp_i <- match(target_data_BART_sub$time, temp_data_BART$startDateTime)
temp_data_BART_sub <- temp_data_BART[BART_temp_i,]

# Specify data (Note: RW = random walk)

data_RW_BART <- list(y = target_data_BART_sub$nee,
                     sd_obs = sd_obs_BART_sub)

data_TempDLM_BART <- list(y = target_data_BART_sub$nee,
                          temp = temp_data_BART_sub$temp,
                          sd_obs = sd_obs_BART_sub)

Specify data from KONZ to be used for model fitting.

# Separate data by site

target_data_KONZ <- filter(target_data, siteID == "KONZ")
temp_data_KONZ <- filter(temp_data, siteID == "KONZ")

# Find data subset w/ fewest number of NAs

min_n_NAs <- subset_length_timesteps  # Initialize scalar for fewest NAs

KONZ_start_i <- NA  # Initialize KONZ starting index

# Loop through data to find data subset with minimum number of NAs for specified subset length

for(i in 1:(length(target_data_KONZ$nee) - subset_length_timesteps + 1)){
  
  num_NAs <- sum(is.na(target_data_KONZ$nee[i:(i + subset_length_timesteps - 1)]))  # Count number of NEE NAs in data subset
  
  if(num_NAs <= min_n_NAs){  # Note: Using <= so that most recent data is used
    
    KONZ_start_i <- i  # Save new starting index
    min_n_NAs <- num_NAs  # Save new lowest number of NAs
    
  }
  
}

# Calculate data subset end index

KONZ_end_i <- KONZ_start_i + subset_length_timesteps - 1

# Subset data

target_data_KONZ_sub <- target_data_KONZ[KONZ_start_i:KONZ_end_i,]

# Interpolate NEE values to estimate sd_obs

nee_interp_KONZ_sub <- approx(x = target_data_KONZ_sub$time[!is.na(target_data_KONZ_sub$nee)],
                              y = target_data_KONZ_sub$nee[!is.na(target_data_KONZ_sub$nee)],
                              xout = target_data_KONZ_sub$time,
                              method = "linear",
                              yleft = mean(target_data_KONZ_sub$nee, na.rm = TRUE),
                              yright = mean(target_data_KONZ_sub$nee, na.rm = TRUE),
                              na.rm = TRUE)

sd_obs_KONZ_sub <- rep(NA, length(target_data_KONZ_sub$nee))  # Initialize sd_obs vector

for(i in 1:length(sd_obs_KONZ_sub)){
  
  if(nee_interp_KONZ_sub$y[i] >= 0){  # Calculate sd_obs if NEE is positive
    
    sd_obs_KONZ_sub[i] <- target_data_KONZ_sub$nee_sd_intercept[1] + target_data_KONZ_sub$ nee_sd_slopeP[1] * nee_interp_KONZ_sub$y[i]  # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
    
  }else{  # Calculate sd_obs if NEE is negative
    
    sd_obs_KONZ_sub[i] <- target_data_KONZ_sub$nee_sd_intercept[1] + target_data_KONZ_sub$ nee_sd_slopeN[1] * nee_interp_KONZ_sub$y[i]
    
  }
  
}

# Match temperature data to subset NEE data

KONZ_temp_i <- match(target_data_KONZ_sub$time, temp_data_KONZ$startDateTime)
temp_data_KONZ_sub <- temp_data_KONZ[KONZ_temp_i,]

# Specify data (Note: RW = random walk)

data_RW_KONZ <- list(y = target_data_KONZ_sub$nee,
                     sd_obs = sd_obs_KONZ_sub)

data_TempDLM_KONZ <- list(y = target_data_KONZ_sub$nee,
                          temp = temp_data_KONZ_sub$temp,
                          sd_obs = sd_obs_KONZ_sub)

Specify data from OSBS to be used for model fitting.

# Separate data by site

target_data_OSBS <- filter(target_data, siteID == "OSBS")
temp_data_OSBS <- filter(temp_data, siteID == "OSBS")

# Remove duplicate rows

target_data_OSBS <- unique(target_data_OSBS)

# Find data subset w/ fewest number of NAs

min_n_NAs <- subset_length_timesteps  # Initialize scalar for fewest NAs

OSBS_start_i <- NA  # Initialize OSBS starting index

# Loop through data to find data subset with minimum number of NAs for specified subset length

for(i in 1:(length(target_data_OSBS$nee) - subset_length_timesteps + 1)){
  
  num_NAs <- sum(is.na(target_data_OSBS$nee[i:(i + subset_length_timesteps - 1)]))  # Count number of NEE NAs in data subset
  
  if(num_NAs <= min_n_NAs){  # Note: Using <= so that most recent data is used
    
    OSBS_start_i <- i  # Save new starting index
    min_n_NAs <- num_NAs  # Save new lowest number of NAs
    
  }
  
}

# Calculate data subset end index

OSBS_end_i <- OSBS_start_i + subset_length_timesteps - 1

# Subset data

target_data_OSBS_sub <- target_data_OSBS[OSBS_start_i:OSBS_end_i,]

# Interpolate NEE values to estimate sd_obs

nee_interp_OSBS_sub <- approx(x = target_data_OSBS_sub$time[!is.na(target_data_OSBS_sub$nee)],
                              y = target_data_OSBS_sub$nee[!is.na(target_data_OSBS_sub$nee)],
                              xout = target_data_OSBS_sub$time,
                              method = "linear",
                              yleft = mean(target_data_OSBS_sub$nee, na.rm = TRUE),
                              yright = mean(target_data_OSBS_sub$nee, na.rm = TRUE),
                              na.rm = TRUE)

sd_obs_OSBS_sub <- rep(NA, length(target_data_OSBS_sub$nee))  # Initialize sd_obs vector

for(i in 1:length(sd_obs_OSBS_sub)){
  
  if(nee_interp_OSBS_sub$y[i] >= 0){  # Calculate sd_obs if NEE is positive
    
    sd_obs_OSBS_sub[i] <- target_data_OSBS_sub$nee_sd_intercept[1] + target_data_OSBS_sub$ nee_sd_slopeP[1] * nee_interp_OSBS_sub$y[i]  # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
    
  }else{  # Calculate sd_obs if NEE is negative
    
    sd_obs_OSBS_sub[i] <- target_data_OSBS_sub$nee_sd_intercept[1] + target_data_OSBS_sub$ nee_sd_slopeN[1] * nee_interp_OSBS_sub$y[i]
    
  }
  
}

# Match temperature data to subset NEE data

OSBS_temp_i <- match(target_data_OSBS_sub$time, temp_data_OSBS$startDateTime)
temp_data_OSBS_sub <- temp_data_OSBS[OSBS_temp_i,]

# Specify data (Note: RW = random walk)

data_RW_OSBS <- list(y = target_data_OSBS_sub$nee,
                     sd_obs = sd_obs_OSBS_sub)

data_TempDLM_OSBS <- list(y = target_data_OSBS_sub$nee,
                          temp = temp_data_OSBS_sub$temp,
                          sd_obs = sd_obs_OSBS_sub)

Specify data from SRER to be used for model fitting.

# Separate data by site

target_data_SRER <- filter(target_data, siteID == "SRER")
temp_data_SRER <- filter(temp_data, siteID == "SRER")

# Remove duplicate rows

target_data_SRER <- unique(target_data_SRER)

# Find data subset w/ fewest number of NAs

min_n_NAs <- subset_length_timesteps  # Initialize scalar for fewest NAs

SRER_start_i <- NA  # Initialize SRER starting index

# Loop through data to find data subset with minimum number of NAs for specified subset length

for(i in 1:(length(target_data_SRER$nee) - subset_length_timesteps + 1)){
  
  num_NAs <- sum(is.na(target_data_SRER$nee[i:(i + subset_length_timesteps - 1)]))  # Count number of NEE NAs in data subset
  
  if(num_NAs <= min_n_NAs){  # Note: Using <= so that most recent data is used
    
    SRER_start_i <- i  # Save new starting index
    min_n_NAs <- num_NAs  # Save new lowest number of NAs
    
  }
  
}

# Calculate data subset end index

SRER_end_i <- SRER_start_i + subset_length_timesteps - 1

# Subset data

target_data_SRER_sub <- target_data_SRER[SRER_start_i:SRER_end_i,]

# Interpolate NEE values to estimate sd_obs

nee_interp_SRER_sub <- approx(x = target_data_SRER_sub$time[!is.na(target_data_SRER_sub$nee)],
                              y = target_data_SRER_sub$nee[!is.na(target_data_SRER_sub$nee)],
                              xout = target_data_SRER_sub$time,
                              method = "linear",
                              yleft = mean(target_data_SRER_sub$nee, na.rm = TRUE),
                              yright = mean(target_data_SRER_sub$nee, na.rm = TRUE),
                              na.rm = TRUE)

sd_obs_SRER_sub <- rep(NA, length(target_data_SRER_sub$nee))  # Initialize sd_obs vector

for(i in 1:length(sd_obs_SRER_sub)){
  
  if(nee_interp_SRER_sub$y[i] >= 0){  # Calculate sd_obs if NEE is positive
    
    sd_obs_SRER_sub[i] <- target_data_SRER_sub$nee_sd_intercept[1] + target_data_SRER_sub$ nee_sd_slopeP[1] * nee_interp_SRER_sub$y[i]  # Note: Since intercept and slope do not change, can just use first value of intercept and slope values
    
  }else{  # Calculate sd_obs if NEE is negative
    
    sd_obs_SRER_sub[i] <- target_data_SRER_sub$nee_sd_intercept[1] + target_data_SRER_sub$ nee_sd_slopeN[1] * nee_interp_SRER_sub$y[i]
    
  }
  
}

# Match temperature data to subset NEE data

SRER_temp_i <- match(target_data_SRER_sub$time, temp_data_SRER$startDateTime)
temp_data_SRER_sub <- temp_data_SRER[SRER_temp_i,]

# Specify data (Note: RW = random walk)

data_RW_SRER <- list(y = target_data_SRER_sub$nee,
                     sd_obs = sd_obs_SRER_sub)

data_TempDLM_SRER <- list(y = target_data_SRER_sub$nee,
                          temp = temp_data_SRER_sub$temp,
                          sd_obs = sd_obs_SRER_sub)

Specify constants for NIMBLE.

# Specify constants

constants_BART <- list(n = length(target_data_BART_sub$nee),
                       x_ic = 0,
                       sd_ic = 1)

constants_KONZ <- list(n = length(target_data_KONZ_sub$nee),
                       x_ic = 0,
                       sd_ic = 1)

constants_OSBS <- list(n = length(target_data_OSBS_sub$nee),
                       x_ic = 0,
                       sd_ic = 1)

constants_SRER <- list(n = length(target_data_SRER_sub$nee),
                       x_ic = 0,
                       sd_ic = 1)

Specify initial conditions for null (random walk) model.

# Specify number of chains

nchains = 3

# Initialize initial condition lists (Note: RW = random walk)

inits_RW_BART <- list()
inits_RW_KONZ <- list()
inits_RW_OSBS <- list()
inits_RW_SRER <- list()

# Generate initial conditions

for(i in 1:nchains){
  
  # BART
  
  y.samp <- sample(nee_interp_BART_sub$y, length(nee_interp_BART_sub$y), replace = TRUE)
  inits_RW_BART[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                             x = nee_interp_BART_sub$y)
  
  # KONZ
  
  y.samp <- sample(nee_interp_KONZ_sub$y, length(nee_interp_KONZ_sub$y), replace = TRUE)
  inits_RW_KONZ[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                             x = nee_interp_KONZ_sub$y)
  
  # OSBS
  
  y.samp <- sample(nee_interp_OSBS_sub$y, length(nee_interp_OSBS_sub$y), replace = TRUE)
  inits_RW_OSBS[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                             x = nee_interp_OSBS_sub$y)
  
  # SRER
  
  y.samp <- sample(nee_interp_SRER_sub$y, length(nee_interp_SRER_sub$y), replace = TRUE)
  inits_RW_SRER[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                             x = nee_interp_SRER_sub$y)
  
}

Specify initial conditions for temperature DLM.

# Initialize initial condition lists

inits_TempDLM_BART <- list()
inits_TempDLM_KONZ <- list()
inits_TempDLM_OSBS <- list()
inits_TempDLM_SRER <- list()

# Generate initial conditions

for(i in 1:nchains){
  
  # BART
  
  y.samp <- sample(nee_interp_BART_sub$y, length(nee_interp_BART_sub$y), replace = TRUE)
  inits_TempDLM_BART[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                                  x = nee_interp_BART_sub$y,
                                  b1 = rnorm(1, mean = 0, sd = 1))
  
  # KONZ
  
  y.samp <- sample(nee_interp_KONZ_sub$y, length(nee_interp_KONZ_sub$y), replace = TRUE)
  inits_TempDLM_KONZ[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                                  x = nee_interp_KONZ_sub$y,
                                  b1 = rnorm(1, mean = 0, sd = 1))
  
  # OSBS
  
  y.samp <- sample(nee_interp_OSBS_sub$y, length(nee_interp_OSBS_sub$y), replace = TRUE)
  inits_TempDLM_OSBS[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                                  x = nee_interp_OSBS_sub$y,
                                  b1 = rnorm(1, mean = 0, sd = 1))
  
  # SRER
  
  y.samp <- sample(nee_interp_SRER_sub$y, length(nee_interp_SRER_sub$y), replace = TRUE)
  inits_TempDLM_SRER[[i]] <- list(sd_add = sd(diff(y.samp), na.rm = TRUE),
                                  x = nee_interp_SRER_sub$y,
                                  b1 = rnorm(1, mean = 0, sd = 1))
  
}

Run NIMBLE for null (random walk) model.

# BART

nimble_out_RW_BART <- nimbleMCMC(code = RandomWalk,
                                 data = data_RW_BART,
                                 inits = inits_RW_BART,
                                 constants = constants_BART,
                                 monitors = c("sd_add",
                                              "x",
                                              "y"),
                                 niter = 11000,
                                 nchains = nchains,
                                 samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
# KONZ

nimble_out_RW_KONZ <- nimbleMCMC(code = RandomWalk,
                                 data = data_RW_KONZ,
                                 inits = inits_RW_KONZ,
                                 constants = constants_KONZ,
                                 monitors = c("sd_add",
                                              "x",
                                              "y"),
                                 niter = 11000,
                                 nchains = nchains,
                                 samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
# OSBS

nimble_out_RW_OSBS <- nimbleMCMC(code = RandomWalk,
                                 data = data_RW_OSBS,
                                 inits = inits_RW_OSBS,
                                 constants = constants_OSBS,
                                 monitors = c("sd_add",
                                              "x",
                                              "y"),
                                 niter = 11000,
                                 nchains = nchains,
                                 samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
# SRER

nimble_out_RW_SRER <- nimbleMCMC(code = RandomWalk,
                                 data = data_RW_SRER,
                                 inits = inits_RW_SRER,
                                 constants = constants_SRER,
                                 monitors = c("sd_add",
                                              "x",
                                              "y"),
                                 niter = 11000,
                                 nchains = nchains,
                                 samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

Run NIMBLE for temperature DLM.

# BART

nimble_out_TempDLM_BART <- nimbleMCMC(code = TempDLM,
                                      data = data_TempDLM_BART,
                                      inits = inits_TempDLM_BART,
                                      constants = constants_BART,
                                      monitors = c("b1",
                                                   "sd_add",
                                                   "x",
                                                   "y"),
                                      niter = 11000,
                                      nchains = nchains,
                                      samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
# KONZ

nimble_out_TempDLM_KONZ <- nimbleMCMC(code = TempDLM,
                                      data = data_TempDLM_KONZ,
                                      inits = inits_TempDLM_KONZ,
                                      constants = constants_KONZ,
                                      monitors = c("b1",
                                                   "sd_add",
                                                   "x",
                                                   "y"),
                                      niter = 11000,
                                      nchains = nchains,
                                      samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
# OSBS

nimble_out_TempDLM_OSBS <- nimbleMCMC(code = TempDLM,
                                      data = data_TempDLM_OSBS,
                                      inits = inits_TempDLM_OSBS,
                                      constants = constants_OSBS,
                                      monitors = c("b1",
                                                   "sd_add",
                                                   "x",
                                                   "y"),
                                      niter = 11000,
                                      nchains = nchains,
                                      samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
# SRER

nimble_out_TempDLM_SRER <- nimbleMCMC(code = TempDLM,
                                      data = data_TempDLM_SRER,
                                      inits = inits_TempDLM_SRER,
                                      constants = constants_SRER,
                                      monitors = c("b1",
                                                   "sd_add",
                                                   "x",
                                                   "y"),
                                      niter = 11000,
                                      nchains = nchains,
                                      samplesAsCodaMCMC = TRUE)
## NAs were detected in model variables: y, logProb_y.
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

Visualize chains and remove burn-in for null (random walk) model.

# BART

plot(nimble_out_RW_BART[, "sd_add"], main = "BART: sd_add")  # Visualize all chains

burnin_RW_BART <- 2000  # Burn-in
nimble_burn_RW_BART <- window(nimble_out_RW_BART, start = burnin_RW_BART)  # Exclude burn-in
plot(nimble_burn_RW_BART[, "sd_add"], main = "BART: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

# KONZ

plot(nimble_out_RW_KONZ[, "sd_add"], main = "KONZ: sd_add")  # Visualize all chains

burnin_RW_KONZ <- 2000  # Burn-in
nimble_burn_RW_KONZ <- window(nimble_out_RW_KONZ, start = burnin_RW_KONZ)  # Exclude burn-in
plot(nimble_burn_RW_KONZ[, "sd_add"], main = "KONZ: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

# OSBS

plot(nimble_out_RW_OSBS[, "sd_add"], main = "OSBS: sd_add")  # Visualize all chains

burnin_RW_OSBS <- 2000  # Burn-in
nimble_burn_RW_OSBS <- window(nimble_out_RW_OSBS, start = burnin_RW_OSBS)  # Exclude burn-in
plot(nimble_burn_RW_OSBS[, "sd_add"], main = "OSBS: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

# SRER

plot(nimble_out_RW_SRER[, "sd_add"], main = "SRER: sd_add")  # Visualize all chains

burnin_RW_SRER <- 2000  # Burn-in
nimble_burn_RW_SRER <- window(nimble_out_RW_SRER, start = burnin_RW_SRER)  # Exclude burn-in
plot(nimble_burn_RW_SRER[, "sd_add"], main = "SRER: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

Visualize chains and remove burn-in for temperature DLM.

# BART

plot(nimble_out_TempDLM_BART[, c("b1")], main = "BART: b1")  # Visualize all chains

plot(nimble_out_TempDLM_BART[, c("sd_add")], main = "BART: sd_add")  # Visualize all chains

burnin_TempDLM_BART <- 2000  # Burn-in
nimble_burn_TempDLM_BART <- window(nimble_out_TempDLM_BART, start = burnin_TempDLM_BART)  # Exclude burn-in
plot(nimble_burn_TempDLM_BART[, c("b1")], main = "BART: b1 (burn-in removed)")  # Visualize chains w/o burn-in

plot(nimble_burn_TempDLM_BART[, c("sd_add")], main = "BART: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

# KONZ

plot(nimble_out_TempDLM_KONZ[, c("b1")], main = "KONZ: b1")  # Visualize all chains

plot(nimble_out_TempDLM_KONZ[, c("sd_add")], main = "KONZ: sd_add")  # Visualize all chains

burnin_TempDLM_KONZ <- 2000  # Burn-in
nimble_burn_TempDLM_KONZ <- window(nimble_out_TempDLM_KONZ, start = burnin_TempDLM_KONZ)  # Exclude burn-in
plot(nimble_burn_TempDLM_KONZ[, c("b1")], main = "KONZ: b1 (burn-in removed)")  # Visualize chains w/o burn-in

plot(nimble_burn_TempDLM_KONZ[, c("sd_add")], main = "KONZ: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

# OSBS

plot(nimble_out_TempDLM_OSBS[, c("b1")], main = "OSBS: b1")  # Visualize all chains

plot(nimble_out_TempDLM_OSBS[, c("sd_add")], main = "OSBS: sd_add")  # Visualize all chains

burnin_TempDLM_OSBS <- 2000  # Burn-in
nimble_burn_TempDLM_OSBS <- window(nimble_out_TempDLM_OSBS, start = burnin_TempDLM_OSBS)  # Exclude burn-in
plot(nimble_burn_TempDLM_OSBS[, c("b1")], main = "OSBS: b1 (burn-in removed)")  # Visualize chains w/o burn-in

plot(nimble_burn_TempDLM_OSBS[, c("sd_add")], main = "OSBS: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

# SRER

plot(nimble_out_TempDLM_SRER[, c("b1")], main = "SRER: b1")  # Visualize all chains

plot(nimble_out_TempDLM_SRER[, c("sd_add")], main = "SRER: sd_add")  # Visualize all chains

burnin_TempDLM_SRER <- 2000  # Burn-in
nimble_burn_TempDLM_SRER <- window(nimble_out_TempDLM_SRER, start = burnin_TempDLM_SRER)  # Exclude burn-in
plot(nimble_burn_TempDLM_SRER[, c("b1")], main = "SRER: b1 (burn-in removed)")  # Visualize chains w/o burn-in

plot(nimble_burn_TempDLM_SRER[, c("sd_add")], main = "SRER: sd_add (burn-in removed)")  # Visualize chains w/o burn-in

Examine diagnostics for null (random walk) model. Since the upper confidence interval for the Brooks-Gelman-Rubin convergence diagnostic is less than 1.1, we can conclude the chains have converged.

# BART

gelman.diag(nimble_burn_RW_BART[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1       1.01
# KONZ

gelman.diag(nimble_burn_RW_KONZ[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1       1.01
# OSBS

gelman.diag(nimble_burn_RW_OSBS[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]       1.01       1.03
# SRER

gelman.diag(nimble_burn_RW_SRER[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]       1.01       1.03

Examine diagnostics for temperature DLM. Since the upper confidence interval for the Brooks-Gelman-Rubin convergence diagnostic is less than 1.1, we can conclude the chains have converged.

# BART

gelman.diag(nimble_burn_TempDLM_BART[, "b1"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]       1.01       1.02
gelman.diag(nimble_burn_TempDLM_BART[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]       1.01       1.03
# KONZ

gelman.diag(nimble_burn_TempDLM_KONZ[, "b1"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1          1
gelman.diag(nimble_burn_TempDLM_KONZ[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1       1.01
# OSBS

gelman.diag(nimble_burn_TempDLM_OSBS[, "b1"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1          1
gelman.diag(nimble_burn_TempDLM_OSBS[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]       1.01       1.04
# SRER

gelman.diag(nimble_burn_TempDLM_SRER[, "b1"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1          1
gelman.diag(nimble_burn_TempDLM_SRER[, "sd_add"])
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## [1,]          1       1.01

Convert null (random walk) model output using tidybayes.

# BART

chain_RW_BART <- nimble_burn_RW_BART %>%
  spread_draws(y[day], x[day], sd_add)

# KONZ

chain_RW_KONZ <- nimble_burn_RW_KONZ %>%
  spread_draws(y[day], x[day], sd_add)

# OSBS

chain_RW_OSBS <- nimble_burn_RW_OSBS %>%
  spread_draws(y[day], x[day], sd_add)

# SRER

chain_RW_SRER <- nimble_burn_RW_SRER %>%
  spread_draws(y[day], x[day], sd_add)

Convert temperature DLM output using tidybayes.

# BART

chain_TempDLM_BART <- nimble_burn_TempDLM_BART %>%
  spread_draws(y[day], x[day], sd_add)

# KONZ

chain_TempDLM_KONZ <- nimble_burn_TempDLM_KONZ %>%
  spread_draws(y[day], x[day], sd_add)

# OSBS

chain_TempDLM_OSBS <- nimble_burn_TempDLM_OSBS %>%
  spread_draws(y[day], x[day], sd_add)

# SRER

chain_TempDLM_SRER <- nimble_burn_TempDLM_SRER %>%
  spread_draws(y[day], x[day], sd_add)

Plot null (random walk) model versus observations.

plot_colors <- c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")

# BART

plot_chain_RW_BART <- chain_RW_BART %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_BART_sub$time)

plot_RW_BART <- ggplot(data = plot_chain_RW_BART, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[1],
              fill = plot_colors[1]) +
  geom_point(data = target_data_BART_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: BART | Time Period: ", substr(target_data_BART_sub$time[1], 1, 10), "to", substr(target_data_BART_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# KONZ

plot_chain_RW_KONZ <- chain_RW_KONZ %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_KONZ_sub$time)

plot_RW_KONZ <- ggplot(data = plot_chain_RW_KONZ, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[2],
              fill = plot_colors[2]) +
  geom_point(data = target_data_KONZ_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: KONZ | Time Period: ", substr(target_data_KONZ_sub$time[1], 1, 10), "to", substr(target_data_KONZ_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# OSBS

plot_chain_RW_OSBS <- chain_RW_OSBS %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_OSBS_sub$time)

plot_RW_OSBS <- ggplot(data = plot_chain_RW_OSBS, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[3],
              fill = plot_colors[3]) +
  geom_point(data = target_data_OSBS_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: OSBS | Time Period: ", substr(target_data_OSBS_sub$time[1], 1, 10), "to", substr(target_data_OSBS_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# SRER

plot_chain_RW_SRER <- chain_RW_SRER %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_SRER_sub$time)

plot_RW_SRER <- ggplot(data = plot_chain_RW_SRER, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[4],
              fill = plot_colors[4]) +
  geom_point(data = target_data_SRER_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: SRER | Time Period: ", substr(target_data_SRER_sub$time[1], 1, 10), "to", substr(target_data_SRER_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

plot_RW_BART

plot_RW_KONZ

plot_RW_OSBS

plot_RW_SRER

Plot temperature DLM versus observations.

# BART

plot_chain_TempDLM_BART <- chain_TempDLM_BART %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_BART_sub$time)

plot_TempDLM_BART <- ggplot(data = plot_chain_TempDLM_BART, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[1],
              fill = plot_colors[1]) +
  geom_point(data = target_data_BART_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: BART | Time Period: ", substr(target_data_BART_sub$time[1], 1, 10), "to", substr(target_data_BART_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# KONZ

plot_chain_TempDLM_KONZ <- chain_TempDLM_KONZ %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_KONZ_sub$time)

plot_TempDLM_KONZ <- ggplot(data = plot_chain_TempDLM_KONZ, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[2],
              fill = plot_colors[2]) +
  geom_point(data = target_data_KONZ_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: KONZ | Time Period: ", substr(target_data_KONZ_sub$time[1], 1, 10), "to", substr(target_data_KONZ_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# OSBS

plot_chain_TempDLM_OSBS <- chain_TempDLM_OSBS %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_OSBS_sub$time)

plot_TempDLM_OSBS <- ggplot(data = plot_chain_TempDLM_OSBS, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[3],
              fill = plot_colors[3]) +
  geom_point(data = target_data_OSBS_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: OSBS | Time Period: ", substr(target_data_OSBS_sub$time[1], 1, 10), "to", substr(target_data_OSBS_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# SRER

plot_chain_TempDLM_SRER <- chain_TempDLM_SRER %>%
  group_by(day) %>%
  summarise(mean = mean(x),
            lower = quantile(x, 0.025),
            upper = quantile(x, 0.975),
            .groups = "drop") %>%
  mutate(time = target_data_SRER_sub$time)

plot_TempDLM_SRER <- ggplot(data = plot_chain_TempDLM_SRER, aes(x = time, y = mean)) +
  geom_line() +
  geom_ribbon(aes(ymin = lower, ymax = upper),
              alpha = 0.2,
              color = plot_colors[4],
              fill = plot_colors[4]) +
  geom_point(data = target_data_SRER_sub, aes(x = time, y = nee), color = "black") +
  labs(title = paste("Site: SRER | Time Period: ", substr(target_data_SRER_sub$time[1], 1, 10), "to", substr(target_data_SRER_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

plot_TempDLM_BART

plot_TempDLM_KONZ

plot_TempDLM_OSBS

plot_TempDLM_SRER

Compare temperature DLM to random walk fit. As can be seen, the two models only slightly differ when there is a relatively large gap in the data.

# BART

comparison_plot_BART <- ggplot() +
  geom_line(data = plot_chain_RW_BART, aes(x = time, y = mean, color = "1")) +
  geom_line(data = plot_chain_TempDLM_BART, aes(x = time, y = mean, color = "2")) +
  geom_ribbon(data = plot_chain_RW_BART, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
  geom_ribbon(data = plot_chain_TempDLM_BART, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
  geom_point(data = target_data_BART_sub, aes(x = time, y = nee), color = "black") +
  scale_fill_manual(name = "Model", values = c("gray", plot_colors[1]), labels = c("Random Walk", "Temp. DLM")) +
  scale_color_manual(name = "Model", values = c("gray", plot_colors[1]), labels = c("Random Walk", "Temp. DLM")) +
  labs(title = paste("Site: BART | Time Period: ", substr(target_data_BART_sub$time[1], 1, 10), "to", substr(target_data_BART_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# KONZ

comparison_plot_KONZ <- ggplot() +
  geom_line(data = plot_chain_RW_KONZ, aes(x = time, y = mean, color = "1")) +
  geom_line(data = plot_chain_TempDLM_KONZ, aes(x = time, y = mean, color = "2")) +
  geom_ribbon(data = plot_chain_RW_KONZ, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
  geom_ribbon(data = plot_chain_TempDLM_KONZ, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
  geom_point(data = target_data_KONZ_sub, aes(x = time, y = nee), color = "black") +
  scale_fill_manual(name = "Model", values = c("gray", plot_colors[2]), labels = c("Random Walk", "Temp. DLM")) +
  scale_color_manual(name = "Model", values = c("gray", plot_colors[2]), labels = c("Random Walk", "Temp. DLM")) +
  labs(title = paste("Site: KONZ | Time Period: ", substr(target_data_KONZ_sub$time[1], 1, 10), "to", substr(target_data_KONZ_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# OSBS

comparison_plot_OSBS <- ggplot() +
  geom_line(data = plot_chain_RW_OSBS, aes(x = time, y = mean, color = "1")) +
  geom_line(data = plot_chain_TempDLM_OSBS, aes(x = time, y = mean, color = "2")) +
  geom_ribbon(data = plot_chain_RW_OSBS, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
  geom_ribbon(data = plot_chain_TempDLM_OSBS, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
  geom_point(data = target_data_OSBS_sub, aes(x = time, y = nee), color = "black") +
  scale_fill_manual(name = "Model", values = c("gray", plot_colors[3]), labels = c("Random Walk", "Temp. DLM")) +
  scale_color_manual(name = "Model", values = c("gray", plot_colors[3]), labels = c("Random Walk", "Temp. DLM")) +
  labs(title = paste("Site: OSBS | Time Period: ", substr(target_data_OSBS_sub$time[1], 1, 10), "to", substr(target_data_OSBS_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

# SRER

comparison_plot_SRER <- ggplot() +
  geom_line(data = plot_chain_RW_SRER, aes(x = time, y = mean, color = "1")) +
  geom_line(data = plot_chain_TempDLM_SRER, aes(x = time, y = mean, color = "2")) +
  geom_ribbon(data = plot_chain_RW_SRER, aes(x = time, ymin = lower, ymax = upper, color = "1", fill = "1"), alpha = 0.2) +
  geom_ribbon(data = plot_chain_TempDLM_SRER, aes(x = time, ymin = lower, ymax = upper, color = "2", fill = "2"), alpha = 0.2) +
  geom_point(data = target_data_SRER_sub, aes(x = time, y = nee), color = "black") +
  scale_fill_manual(name = "Model", values = c("gray", plot_colors[4]), labels = c("Random Walk", "Temp. DLM")) +
  scale_color_manual(name = "Model", values = c("gray", plot_colors[4]), labels = c("Random Walk", "Temp. DLM")) +
  labs(title = paste("Site: SRER | Time Period: ", substr(target_data_SRER_sub$time[1], 1, 10), "to", substr(target_data_SRER_sub$time[subset_length_timesteps], 1, 10)),
       x = "Time (UTC)",
       y = expression(paste("NEE (",mu, "mol" ~ CO[2] ~ m^-2 ~ s^-1,")"))) +
  theme_bw()

comparison_plot_BART

comparison_plot_KONZ

comparison_plot_OSBS

comparison_plot_SRER